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Abstract 

■ Excited states of Bose-Einstein condensates are considered in the semi-classical (Thomas- 

Qh ' Fermi) limit of the Gross-Pitaevskii equation with repulsive inter-atomic interactions and a 

harmonic potential. The relative dynamics of dark solitons (density dips on the localized 
condensate) with respect to the harmonic potential and to each other is approximated using 
the averaged Lagrangian method. This permits a complete characterization of the equilibrium 
positions of the dark solitons as a function of the chemical potential parameter. It also yields 
an analytical handle on the oscillation frequencies of dark solitons around such equilibria. 
^ ' The asymptotic predictions are generalized for an arbitrary number of dark solitons and are 

, corroborated by numerical computations for 2- and 3-soliton configurations. 

(N 

^ ■ 1 Introduction 

O _ 

The defocusing nonlinear Schrodinger equation is a prototypical model for a variety of different 
settings including nonlinear optics, liquids, mechanical systems, and magnetic films, among others. 
In one spatial dimension, its prototypical excitation is the dark soliton, i.e., a localized density 
dip on a continuous- wave background (carrying also a phase jump). 

One of the major areas where the description of dark solitons with a mean-field model (also 
known as the Gross-Pitaevskii equation) has been the physics of atomic Bose-Einstein condensates 
(BECs) |13|. I14j. There, the repulsive inter-atomic interactions can be accurately captured by an 
effective nonlinear self-action [4J. A considerable volume of experimental work has conclusively 
demonstrated the relevance of such nonlinear waveforms within harmonically confined condensates. 
Although in earlier works, such coherent structures were dynamically or thermally unstable [21 
[3], more recent work has overcome such limitations [T5| [TfJl I17j . This has been achieved by 
working at sufficiently low temperatures (of the order of lOnK) and for strongly confined in the 
transverse directions, cigar-shaped BECs. Furthermore, in these recent experiments, the nature 
of the generation process (e.g., by interference of two independent BECs |1 5| . 117 } [T8]. or through 
interaction of the BEC with an appropriate light pulse [E]), it has been possible to produce two 
or more dark solitons on the background of a localized condensate. In principle, the resulting 
number of dark solitons can be chosen at will, as indicated in [17j . 

These recent developments prompt us to examine the dynamics of dark solitons which are har- 
monically confined within localized repulsive Bose-Einstein condensates. These can be thought of 
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as density dips that arise in nonlinear variants of the excited states of the quantum harmonic os- 
cillator pQ. The study of the equilibrium positions and near-equilibrium dynamics of these density 
dips is the principal theme of the present contribution. In particular, using a Lagrangian (vari- 
ational) approach, we compute the asymptotic dependence on the chemical potential parameter 
both for equilibrium positions of dark solitons and for their oscillation frequencies around such 
equilibria. 

This article is organized as follows. In Section 2, we present the general mathematical setup 
of the problem. Section 3 examines the single soliton case, Section 4 extends considerations to 
2-solitons, and Section 5 generalizes the results to an arbitrary number of m-solitons for m ^ 2. 
Section 6 compares our asymptotic predictions to numerical computations and suggests some 
interesting directions for further study. 



2 Mathematical Setup 

Let us start with the Gross-Pitaevskii equation with a harmonic potential and repulsive nonlinear 
interactions ^ ^ 

iv T = --v^ + -fv+\v\ 2 v - fiv, (1) 

where i>(£,r) : M x 8 -> C is the wave function and /i G M. represents the chemical potential 
(and is physically associated with the number of atoms in the condensate). We are interested in 
localized modes of the Gross-Pitaevskii equation in the limit /i — > oo, which is associated with the 
semi-classical or Thomas-Fermi limit. Using the scaling transformation, 

v(£,t) = /// 2 u(x,i), C = (2 / u) 1 / 2 x, r = 2t, (2) 

the Gross-Pitaevskii equation ([1]) is transformed to the semi-classical form 

ieut + £ 2 u xx + (1 — x 2 — \u\ 2 )u = 0, (3) 

where u(x, t) : K x R — > C is a new wave function and e = (2/j 1 )~ 1 is a small parameter. 
Let r] £ be a real positive solution of the stationary problem 

s 2 r]' l !(x) + (l-x 2 -T] 2 (x))r] £ (x) = 0, x G R. (4) 

Main results of Ignat & Millot [8l [9] and Gallo & Pelinovsky [6] state that for any sufficiently 
small e > there exists a smooth solution rj e G C°°(R) that decays to zero as \x\ — > oo faster 
than any exponential function. The ground state converges pointwise as e — ► to the compact 
Thomas-Fermi cloud 

, , , > / (1 - x 2 ) 1 / 2 , for \x\ < 1, , , 

r ?0 (x):=hm ??£ (x) = | ^ for \x\ > 1. (5) 

Useful properties of the ground state r] £ for sufficiently small e > are summarized as follows: 
• For any compact subset K G (—1,1), there is Ck > such that 

he ~ VoWcHK) < C K e 2 . (6) 
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• There is C > such that 

||% - 7?o||loc < Ce 1/3 , II^Hlcc < Ce" 1 / 3 , \\r£\\ L ~, «S Ce' 1 . (7) 

We shall consider excited states of the Gross-Pitaevskii equation ([3]), which are non-positive 
solutions of the stationary problem 

e 2 u"(x) + (l-x 2 -u 2 £ (x))u £ (x) = 0, i£R. (8) 

The excited states can be classified by the number m of zeros of u £ (x) on R. A unique solution 
with m zeros exists near e = e m by the local bifurcation theory [12j . where e m = 1+ 1 2m , m G N. 
Because of the symmetry of the harmonic potential, the m-th excited state u £ (x) is even on R 
for even m € N and odd on R for odd m £ N. The m-th excited state is continued for e < e m 
numerically by Zezyulin et al. |19j . 

In our work we shall apply variational approximations [TT] to study relative dynamics of dark 
solitons (localized solutions with nonzero boundary conditions on the background of the positive 
ground state rj E ) with respect to the harmonic potential and to each other. In particular, we 
obtain results on existence and spectral stability of the excited states from analysis of equilibrium 
positions of dark solitons and their oscillation frequencies near such equilibrium. To enable this 
formalism, we substitute 

u(x,t) = r] E (x)v(x, t) 
to the Gross-Pitaevskii equation (|3|) and find an equivalent equation 

ierj 2 e v t + e 2 (rj 2 v x ) x + r£(l - \v\ 2 )v = 0. (9) 

Excited states are solutions of the stationary equation 

S2 TX (^)^)) + " V m^Wm(x) =0, X <E R, (10) 

which have exactly m zeros on R and satisfy the boundary conditions 

lim V m (x) = (±l) m , m £ N. 

Solutions of the stationary Gross-Pitaevskii equation (|10j) are critical points of the energy func- 
tional 

A( V ) = e 2 [ Ve(x)K\ 2 dx + \l 4{x){l - \v\ 2 ) 2 dx. (11) 

in the sense of 4=-| v= y m = 0. The time-dependent Gross-Pitaevskii equation ((H) follows from the 
Lagrangian function L(y) = K(v) + A(v), where 
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K(v) = -e / r] £ (x)(vv t - vv t )dx, (12) 



by means of the Euler-Lagrange equations 

5L d 5L _ 
5v dt 5vt 

In what follows, we obtain variational approximations for time-dependent solutions near the ex- 
cited states V m {x) for m = 1, m = 2, and in the general case m ^ 2. We also compare these 
approximations with numerical results for m = 2 and m = 3. 
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3 1-soliton (m = 1) 

Let us consider the dark soliton 

vi(x,t) = A(t) tanh(e' 1 B(t)(x - a(t))) +ib(t), A > 0, B > 0, a £ M, b e R, (13) 

as an ansatz for the Lagrangian L(y). The motivation for this choice originates from the fact that 
(fl~3l) is an exact solution of (jHJ) if % = 1 under constraints 

^4 = ^/1- 6 2 , 5 = ^Vl-& 2 , a = a + \/26t, b = b , 
V2 

where ao £ R and &o £ ( — 1) 1) are arbitrary t-independent parameters. In view of the relation 

\ Vl \ 2 = A 2 + b 2 - A 2 sech 2 (e^Btyix - a(t))) , 

it is clear that a is a center of the dark soliton, b its speed, A determines its amplitude, and B 
determines its width. If the dark soliton is placed inside the confinement of the compact Thomas- 
Fermi cloud ©, then the constraint a G (—1,1) has to be added. 

When r/ £ / 1, the trial function (I13p is no longer an exact solution of ([9]) but it becomes the best 
approximate solution if parameters (A, B, a, b) are chosen from the Euler-Lagrange equations of the 
averaged Lagrangian L\(A, B, a, b) = L(v\). This variational method provides a useful qualitative 
approximation to physicists for understanding the dynamics of dark solitons under perturbations 
|llj . Unlike the work of |11| . we do not need to renormalize the Lagrangian function L(v) thanks 
to the rapidly decaying weight function r] 2 (x) under the integration sign in (|lip - (|12l) . 



Let us choose A = yl — b 2 to satisfy the boundary conditions 

lim \vi(x,t)\ = 1 for all tel. 

x— »±oo 

Substitution of ansatz (]13p to L(v) and integration in R results in the effective Lagrangian 

L(yx) = £ ^ f r] 2 (x) tanh(z)dx + b\/ 1 — b 2 Bd f 7? 2 (x)sech 2 (z)dx 
V 1 - b 2 Jr Jr 

-eby/l - b 2 BB~ l / rj 2 (x)zsech 2 (z)dx + (1 - b 2 )B 2 I T/ 2 (x)sech 4 (z)dx 
Jr Jr 

+W~b 2 ) 2 [ vt(x)sedi A (z)dx, (14) 
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where z = e~ 1 B(x — a). Note the pointwise limits 

lim tanh(z) = sign(x — a), lim sech 2 (z) = 0, x G R\{0}, (15) 

which show that lim e ^Q L(vi) = 0. The value of L{v\) in the limit of e — > is computed in the 
following lemma. 

Lemma 1 Assume that B > and a £ (—1, 1). Then, 

Li:=lim ^M = -- -j— (q - -a 3 ) + bVT^jl - a 2 )d 



e- 



+ o 2e VT^ti 2 ' 3 



+^(l-a 2 )(l-6 2 )5 + ^(l-a 2 ) 2 (l-6 2 ) 2 . 
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Proof. Thanks to the limit ([5]), the pointwise bound (|15p . and the Dominated Convergence 
Theorem, we have 



lim / i]e( x ) ta.vh(z)dx = [ (1 — x 2 )sign(x — a)dx = —2a H — c 
£ ^°Jr J -i 3 

To compute the remaining four integrals in (|14p . we use the change of variables x — > z, so that 

e -1 ^ / ^(a;)sech 2 (2;)dx = / rj^ (a + ezB^ 1 ) sech 2 (z)dz 
Jr Jr 

= / rjQ (a + ezB^ 1 ) sech 2 (z)dz + e 1/3 R £ , B ,a{ z ) sech2 ( z )dz, 

J z_ Jr 

where z± = e~ 1 B(±l — a) and the reminder term satisfies the uniform bound H-Re.^alU 00 ^ C for 
some C > 0, thanks to the first bound ([7]). As a result, the second term does not contribute to 
the limit e — > 0. To deal with the first term, we decompose the integral into three parts 

fz+ rz+ pz + 

{I -a 2 ) sech 2 (z)dz -2ECIB- 1 zsech 2 (z)dz - e 2 B~ 2 z 2 sech 2 (z)dz. 

J Z— J Z— J Z- 

We recall that the integral 

\ z k sech 2 (z)dz, k ^ 0, 

Jae- 1 

is exponentially small in e if a > is e-independent. As a result, the second and third terms do 
not contribute to the limit e —* 0, while the first term gives 



lime 1 B / r] 2 (x)sech 2 (z)dx = (1 — a 2 ) / sech' 
Jr Jr 



(z)dz = 2(1 - a 2 ). 



The remaining three integrals in (|14p are computed similarly to the second integral in (|14p and 
give 

lime -1 !? / i] 2 (x)zsech. 2 (z)dx = (1 — a 2 ) / zsech 2 (z)dz = 0, 
Jr Jr 

lime -1 !? / r] 2 (x)sech A (z)dx = (1 — a 2 ) J sech A (z)dz = -(1 — a 2 ), 
Jr Jr 3 

lime^B / r] 4 (x)sedi A (z)dx = (1 - a 2 ) 2 / sech 4 (z)dz = -(1 - a 2 ) 2 . 
Jr Jr 3 

Combining all individual computations gives the result for L\. ■ 

Since B is absent in L\ := L\(a,b, B), variation of L\ with respect to B gives an algebraic 
equation on B with the exact solution 

B= -Wl-aVl-fe 2 - 
Eliminating from L\(a, b, B), we simplify the effective Lagrangian to the form 



Ll (a, b) = ^(1 - a 2 ) 3 / 2 (l - 6 2 ) 3 / 2 - 2y^V 2 b(a - \a*) + | 



(a - ^a 3 )bVl ~ b 2 



where the last term is the full derivative. Since adding a full derivative does not change the Euler- 
Lagrange equations, the last term can be dropped from L\. Variation with respect to a and b give 
the following system of equations 

r I o ; V2a(l-b 2 ) 

a = V2Vl-a 2 b, b = - -^ = 1 , 

Vl-o 2 

which is equivalent to the linear oscillator equation 

a + 2a = 0. 

The critical point (a, b) = (0,0) corresponds to the solution V\ of the stationary equation (fTUj) , 
Oscillations near the critical point with frequency v2 corresponds to the oscillations of the dark 
soliton V\ relative to the positive ground state rj e in the Thomas-Fermi limit e — > 0; see e.g. |10] 
and references therein. This frequency was found to be the smallest nonzero frequency in the 
spectrum of the spectral stability problem associated with the first excited state, see Fig. 2 in 

4 2-solitons (m = 2) 

Let us now consider a superposition of two dark solitons 

v 2 (x,t) = tanh(e~ 1 5i(t)(x-ai(t))) 

x [A 2 (t) tanh (e- 1 B 2 (t)(x - a 2 (t))) + ib 2 {t)] , (16) 

where we shall use the relations for the individual dark solitons 

A 3 = ^l £. = _L0T^IT^ i = 1]2 . 

In-phase oscillations of two dark solitons are very similar to the oscillations of one dark soliton and 
have the same frequency, as we will show in Section [5j Therefore, we shall consider out-of-phase 
oscillations of two dark solitons and choose 

a% = —a, a 2 = a, bi = —b, b 2 = b, 

with a S (0, 1) and b £ R. Substitution of v 2 to A(v) gives 

A( V2 ) = A 2 B 2 [ rilix) [sech 4 (z+) +sech 4 (z_) - 26 2 sech 2 (z+)sech 2 (z_) 

— ^4 2 sech 2 (2; + )sech 2 (2;_) (sech 2 (z + ) +sech 2 (z„) — 2 tanh(z + ) tanh(z„))] dx 

+ \a a ( 7]^(x) [sech 4 (z + ) + sech 4 (z_) + 2sech 2 (z + )sech 2 (z_) 

-2A 2 sech 2 (z+)sech 2 (z_) (sech 2 (z + ) +sech 2 (z_)) + A 4 sech 4 (z + )sech 4 (z_)] dx, 

where z± = e~ 1 B{x ± a). The integrals that only depend on z + or z_ are computed similarly 
to the case of 1-soliton. The overlapping integrals that depend on both z + and Z- are computed 
under the apriori assumption 

a^de 1/6 , e - iBae ^ < C 2 e 2 log(e), (17) 



6 



for some C\,C2 > and sufficiently small e > 0. As we will see later, the apriori assumption 
allows us to recover the equilibrium state of two dark solitons and to study perturbations near the 
equilibrium. 

After simplifications, one can write 

A M V 2) A A A 

A2 := ^ — A + + A - + A over i a p, 

where 

A 2 B 2 f o , , . „ , , . A 4 



A± :-- 



2s 
and 

12 D 2 



/7/ 2 (x)sech 4 (z-|-)(ix + — f f] A (x)sech 4: (z±)da 



Aoverlap = ~ — / ?7 e (x)sech 2 (z + )sech 2 (z- ) 

x [2b 2 + yl 2 (sech 2 (z+) + sech 2 (z_) - 2tanh(z + ) tanh(z_))] dx 
+— / ^(x)sech 2 (z_|_)sech 2 (2;_) 

x [2 - 2A 2 (sech 2 (z+) + sech 2 (z_)) + ,4 4 sech 2 (z + )sech 2 (z_)] dx. 



The terms A± are the potential energies of the individual dark solitons and the term A over i ap 
contains overlapping integrals. By Lemma [H we have 

A± = 4 (1 - y(l-^ +0(el/3) 
3v 2 

The overlapping integrals for small e are computed in the following lemma. 
Lemma 2 Assume that a £ (0, 1) satisfies b € K, and 

A = Vl-6 2 , £ = ^\A-a 2 \A-fr 2 - 
V2 

A over i ap = -8V2(1 - a 2 ) 3 / 2 (l - 6 2 ) 5 / 2 e-^ 1 ( 1 + O^ 3 ;' 



Proof. To compute the overlapping integrals, we use the symmetry of the integrand and the 
change of variables x — ► Z- . The first overlapping integral in A over i ap is given by 

s~ 1 B / n 2 (a;)sech 2 (z + )sech 2 (z_)(ix = 2 / rj 2 (a + ezB' 1 ) sech 2 (z)sech 2 (z + 2Bae~ 1 )dz, 

where z = z~ . Similarly to the proof of Lemma [H we break the integral into four parts 

rB{\-a)e- x 

2(1 -a 2 ) / sech 2 (z)sech 2 (z + 2 J Bae^ 1 )dz 

J-Bae- 1 
rB{X-a)e~ x 

-4asB~ L / zsech 2 (z)sech 2 (z + 2Bae^ 1 )dz 

J-Bae- 1 
rB(l-a)e- 1 

-2e 2 B- 2 / z 2 sech 2 (z)sech 2 (z + 2 J Bae" 1 )(iz 

J-Bae- 1 

/oo 
R £ ,B,a( z ) sech2 ( z ) sech2 ( z + 2Bas~ L )dz, 
-Bae- 1 
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where the reminder term satisfies the bound [|i?E s a||i°° ^ C f° r some C > 0, thanks to the bound 
(J7j) . The first part gives the leading order of the integral according to the explicit calculation 



/vB(l-a)e -1 

Ji = / sech 2 (z)sech 2 (z + 2Bae~ 1 )dz 

J-Bae- 1 



Bae- 1 r B(l-a)e- 1 \ ^iz-iBae' 1 



We have 
and 
so that 



' ' -Baz-^JBae-i J (1 + e" 22 ) 2 (l + e^"^ 1 )2 ^ 



< e -^- 4Ba£ 1 < e~ 2BaE \ z ^ -Bae- 1 , 

e -B(l-a)e-^ e -Bae-i fl ^ ^1/6 



= Se" 4 ^ 1 (2Bae- 1 - l) (l + O (e" 2 ^' 1 
The second part of the overlapping integral is computed from the explicit computation 

f .B(l-a)e- 1 



/±S(i— a)e 
zsech 2 (z)sech 2 (z + 2Bae~ 1 )dz 
-Bae- 1 

/ r-Bae- 1 rB{X-a)e- x \ 

ae / + / zsech 2 (z)sech 2 (z + 2Bae~ 1 )dz 

\J-Bae- 1 J Bae- 1 ) 



= 0{a 2 h) + O [e- &Bae J = 0(a 2 h), 
The last two parts of the overlapping integrals are computed similarly and yield 

i-B(l-a)e- 1 

h = e 2 / z 2 sech 2 (z)sech 2 (z + 2Bae- 1 )dz = 0(a 2 I 1 ), 

J-Bae- 1 



h = e 1/3 



/oo 
-R e ,s,a(2)sech 2 (z)sech 2 (z + 2Bae~ 1 )dz = 0(e 1/3 h). 
Bae- 1 

Under the assumption (|17|). we have 

e- 2Bae ^ =0{e\og 1,2 {e)) and a 2 = 0(e 1/3 ), 

so that we finally obtain 

e- x B f 7 1 2 {x)sech 2 (z + )sech 2 (z_)dx = 16(1 - a 2 )e~^ Ba£ ' 1 (2Bae~ 1 - l) f 1 + 0(e 1/3 ) 
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Similarly, we compute the other overlapping integrals in A over i ap as follows: 

s~ l B / ?7^(x)sech 2 (z_|_)sech 2 (z_) (sech 2 (z + ) + sech 2 (z_)) dx 
Jr 

= ^(l-a 2 )e- 4Ba£ " 1 (l + OCe 1 / 8 )), 
e~ x B / ry 2 (x)sech 2 (z+)sech 2 (z_)tanh(z + )tanh(2;_)cZa; 

Jm. 

= 32(1 - a 2 )e- 4Sae_1 {-Bae~ l + l) (l + ©(e 1 / 3 )) , 

and 

e^B J r ? 4 (x)sech 4 (z + )sech 4 (z_) ( ix = 512(1 -a 2 ) 2 e- 8Ba£_1 ^ae" 1 - (l + 0(^/3- 

Combining these computations together, we obtain the expression for A over i ap . ■ 

Variations of A^(a, b) define critical points that correspond to the solution ^(z) of the sta- 
tionary equation (jlOp . Since A2 is even in b 6 R, the set of critical points includes 6 = 0. Note 
that V2(x,t) in (|16|) is real if b = 0, which agree with V^(x) being real-valued. 

Since A + + A_ is even in a and the overlapping integral is small under assumption ([T7]) . 
variation of A2(a,0) in a gives a root finding problem 

- 4V2ea (l + 0(e 1/3 )) + 32e~ 2v/2ae ~ 1 (l + ©(e 1 / 3 )) = 0. (18) 



The asymptotic analysis of the roots of the nonlinear equation (|18p in the following lemma shows 
that the apriori assumption (|17p is indeed satisfied. 



Lemma 3 For sufficiently small e > 0, there exists a simple root of the nonlinear equation j!8\ ) 
in the neighborhood of 0, which is expanded by 

o = -^^-log(e)-ilog|log(e)| + |log(2) + o(l)) as e -» 0. (19) 



Proof. Taking a natural logarithm of the nonlinear equation (|18p , we obtain 

2v / 2a + elog(a) = -elog(e) + ^elog(2) + 0(e 4/3 ). 
Let a = — -j=elog(e)U and rewrite the problem for U: 

_ \og(U) log|log(e)| _ 31og(2) / 1/3 

21og(e) + 21og(e) 21og(e) V + U{£ ' 



(20) 



By the Implicit Function Theorem applied to equation (|20p . existence of a unique root U{e) in a 
one-sided neighborhood of e > is proved, where U(e) is continuous in e > and rirn e io U(e) = 1. 
To estimate the remainder term for \U(e) — 1|, one can further decompose 

p _ 1 + ispsgL 

21og(e) 
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and rewrite the problem for V: 



V 



lQR(l + ^gf(l + F) 

log | log (e) | 



31og(2) 
log | log(e) 



1 + 0(e 1/3 )) . 



(21) 



By the Implicit Function Theorem applied again to equation (|2ip . existence of a unique root 
V(e) in a one-sided neighborhood of e > is proved, where V(e) is continuous in e > and 
lim e jo = 0. Substitution of U back to formula for a gives (fT9j) , ■ 

By Lemma O we can study temporal dynamics of two dark solitons near the bound state that 
corresponds to a small root of the nonlinear equation (|18p . 

To proceed with time-derivative terms, we substitute (fT6|) to the kinetic part K(v) in (fT2|) and 
find that 

— A+ + + -^overlap, 



#2 := 



2t 



where 



1T± 



and 



A', 



2^1^62 

bVT^B 
' 25 



1 ■ 



r] 2 (x) tanh(z_|-)da; + 



2s 



rj 2 (x)sech 2 (z± ) dx 



rj 2 (x)z± sech 2 {z±)dx 



overlap 



6(1 — 6 2 ) 1/2 / if £ (x) (tanh(z + )sech 2 (z_) - tanh(z_)sech 2 (z+)) dx 
2 Jr 

-e~ 1 b(l-b 2 ) 3/2 (Ba + Ba) [ r ] 2 (x)sech 2 (z + )sech 2 (z-)dx. 

Jr 

The terms K± are the kinetic energies of the individual dark solitons and the term -Koverlap contains 
overlapping integrals. By Lemma [H we have 



1 



d 



lim( J fC+ + K-) = -4\A - b 2 b{a - -a*) + 2 , 
e^O 3 dt 



(a - \a z )b^l - b 2 



The overlapping integrals for small e are estimated in the following lemma. 



Lemma 4 Assume that a G (0, 1) satisfies (11), b G R, and 



v2 



Then, 



^overlap = 2e6(l - b^B' 1 (1 - G 2 ) (l + O (e 1 / 3 

-166(1 - 6 2 ) 3 / 2 (d + B^Ba){l - a^e" 45 ^" 1 (2Bae~ 1 - l) (l + 0(e 1/3 )) . 



Proof. The first and second terms in i^overiap are estimated similarly to the proof of Lemma [2j 
Note that the first term disappears in the limit e — » 0. ■ 
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To obtain effective dynamical equations on (a, b) valid in the domain specified by assumption 
()17p . we expand L 2 {a, b) = in the quadratic form in (a, b) and apply the limit e —* to all 
but the overlapping integrals. As a result, the reduced effective Lagrangian L2(a, b) takes the form 



L 2 (a,b) ~ "I"! 1 2 ' 



+ a 2 ) + 0(6 2 + a 2 ) 2 ) - 4ab (l + 0(6 2 + a 2 )) 
-S^ 2 ^ 1 ^ ^ 2 )) (1 + 0(6 2 + a 2 )) . 
In variables (a,b), the Euler-Lagrange equations at the leading order become 

a = V2b, b = -V2a + 8£- 1 e- 2 ^ as ~\ 
or, equivalently, recover the nonlinear oscillator equation 



t— i 2\/2a 

a + 2a = 8V2e e ~. 



The equilibrium state is given by the root ao(e) of the nonlinear equation (|18p . This equilibrium 
state is a center and linear oscillations near the center satisfy 

5 + ujq5 = 0, 

where 5 = a — ao(e) and 

^(e) = 2 + SL-aVSooMe-i = 2 + l^f) 

= -41og(e)-21og|log(e)|+2 + 61og(2) + o(l), as e -» 0, (22) 



thanks to Lemma [3l We note that the frequency ujq(e) of out-of-phase oscillations of two dark 
solitons grows in the limit e — > 0. This property will be further discussed in Section El 

5 m-solitons with m ^ 2 

We extrapolate the results of the previous section to the case of m-solitons with m ^ 2. The 
general superposition of m dark solitons is substituted in the form 

m 

) = Yl (Aj(t) tanh (e^Bji^ix - aj (t))) + ibj(t)) , (23) 
j=l 

where 

1 



A j = ^l-b], B 3 = —^l-a^l-b], j€{l,2,...,m} 
Under the same assumptions of 

\aj\^Ce 1/6 , j E {1, 2, m} 

and 

g-^+x-a^e- 1 ^C£ 2 log(e), je{l,2,..,m-l}, 
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for some C > 0, we reduce the effective Lagrangian L m := to the leading order 

m m m— 1 

j=i J= i j=i 

where only the quadratic terms in (a,j,bj) and only the pairwise interaction potentials are taken 
into account. Using the Euler-Lagrange equations, we obtain 

aj = y/2bj, bj = -V2aj - ^-V^fe+i-a,)^ 1 _ g-v^-a^)^^ ^ ■ g ^ (34) 

where boundary conditions ao = — 00 and a m+ \ = 00 must be used. The center of mass (a) = 
— 5^j=i a j satisfies the linear oscillator equation 

(a) + 2(a) = 0, (25) 

which recovers the frequency of oscillations of a 1 dark soliton in Section [3J Let us introduce the 
set of normal coordinates 

xj = V2(a j+ i -a,j)e~ l , j G {1,2, ...,m - 1}, 

and rewrite system (|24p in the scalar form 

xj + 2xj + I6e~ 2 (e~ x ^ +1 - 2e - ^ ' + e'^- 1 ) = 0, j £ {1, 2, m — 1}, (26) 

where the boundary conditions are now xq = x m = 00. System (I26D is known as the Toda lattice 
with nonzero masses, which is not integrable by inverse scattering (unlike its counterpart with 
zero masses). We are only interested in existence of critical points in the Toda lattice and in the 
distribution of eigenvalues in the linearization around the critical points. 

Critical points of the Toda lattice (I26p are defined by solutions of system of algebraic equations 

2xj + 16e~ 2 (e~^'+ 1 - 2e - ^ + e~ x ^ 1 ) = 0, j £ {1, 2, m — 1}. (27) 

Let the (m — 1) x (m — 1) matrix A be given by 



2 


-1 





••• 








-1 


2 


-1 













-1 


2 


_1 ... 

















... 


-1 


2 



Matrix A arises in the central-difference approximation of the second derivative subject to the 
Dirichlet boundary conditions. It is strictly positive and thus invertible. The system of algebraic 
equations (|27p can be written in the matrix-vector form 

2 2 

Ae~ x = — x e~ x = — A~ x x. (28) 

8 8 

Solutions of system (|28p in the limit e — > are analyzed in the following lemma. 
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Lemma 5 For sufficiently small e > 0, there exists a unique solution of system \28\) in the 
neighborhood of oo, which is expanded by 

x = -21og(e)l-log|log(e)|l + 21og(2)l-log(A" 1 l)+o(l), as e 0, (29) 

where 1 = [1, 1, 1] T G R m_1 . 

Proof. Applying the natural logarithm to system (I28p . we rewrite the system as follows 

x = -21og(e)l + 31og(2)l - log (A _1 x) . 

Repeating the proof of Lemma El we find the desired expansion (i29j) . ■ 

Back to the physical variables (oi, a m ), the result of Lemma[5]implies that the coordinates of 
dark solitons are centered (a) = and distributed with nearly equal spacing as e — > 0. Linearizing 
the Toda lattice (|26|) about the root of system (|27|) . we obtain the linear eigenvalue problem 

(2 - u/% - I6e~ 2 (e-^'+^+i - 2e~ x ^ 3 + e'^-^j-i) = 0, j G {1, 2, m - 1}, (30) 

where £o and £ m are not determined because the coefficients in front of £o and £ m are zero. Using 
the representation (|28p . we rewrite the linear eigenvalue problem in the form 

(2_ w %-2((A- 1 x) i+1 ^ +1 -2(A- 1 x) J ^ + (A- 1 x)^ 1 ^_ 1 ) =0, j G {1, 2, m — 1}. (31) 

Frequencies of oscillations are analyzed in the limit e — ► in the following lemma. 

Lemma 6 For sufficiently small e > 0, (m— 1) eigenvalues of the linear problem Ii31\) are expanded 
by 

to 2 = 2 + (-41og(e) - 2 log | log(e)| + 41og(2)) Q 2 + 0(1), (32) 
where Cl 2 G |l, 3, 6, lMl?L_U. j ari j m ^ 2. 

Proof. Let SI 2 be eigenvalues of the reduced eigenvalue problem 

Sl 2 £j + Vj + i€ j+1 -2vj€ j + Vj- 1 £j-i = 0, j G {1, 2, m — 1}, (33) 

where v = A _1 l G JR m_1 . We will show that all eigenvalues of the reduced eigenvalue problem 
(f33|) are simple and given explicitly by fi 2 G {l,3,6,.., nfcil}. If this is the case, the asymptotic 
expansion (j29[) and the regular perturbation theory for the matrix eigenvalue problem (|31f) imply 
that 

\oj 2 - 2 + (41og(e) + 21og|log(e)| -41og(2))^ 2 | = O(l), as e 0, 
for each eigenvalue Q 2 . 

To obtain the exact distribution of eigenvalues of the reduced eigenvalue problem ([33]) , we will 
find the vector v explicitly. The components of v satisfy the Dirichlet problem for second-order 
difference equations 

2vj - v j+ i - Vj-i = 1, j G {1, 2, m - 1}, 
subject to vq = v m = 0. The exact solution of this problem is 

Vj = ^j{ m -j), j G {1,2, ...,m- 1}. 
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Let k = j - so that k G X m := {-?r + 1, -f + 2, f - 1}. Note that J TO includes integer 
values for even m and half-integer values for odd m. Denote Ck = £i, an d A = 2f2 2 and rewrite the 
reduced eigenvalue problem (|33|) in the following explicit form 

ACfc = - k 2 ^j {2( k - (k+i ~ Ck-i) + 2k (Cfe+i - Ck-i) + (Cfc+i + Cfc-i) , * e X m - (34) 

First, we consider the problem (134p for all /c G Z with a fixed m ^ 2 and prove that there exists 
a basis of eigenvectors £ G {P n }neN i n the space of analytic functions on Z for an infinite set 
of eigenvalues A G {(n + l)(n + 2)} ng N , where No := {0, 1, 2, ...}. The corresponding eigenvector 
£ = P n for each eigenvalue A = (n + l)(n + 2) is given by the polynomial P n (k) in the form 

Q k = P n (k):=k n + c 1 k n - 1 + c 2 k n - 2 + ... + c n , k£Z, (35) 

with uniquely determined coefficients {c\, C2, c n ). To show this, we note that if C £ Pni where 
V n is the vector space of polynomials of degree n, then the vector field of the eigenvalue problem 
(I34p belongs to V n . This follows from the fact that if £ G "P n , then 

(2Cfc-Cfc+i-C*-i)ePn-2, (Cfc+i - Ck-i) e TVi, (Cfc+i + Cfc-i) e Pn- (36) 

Substituting the representation ()35[) to the linear eigenvalue problem (|34p . we collect coefficients 
in front of k n to find that A = (n + l)(n + 2) and the coefficients in front of /c n_1 , k n ~ 2 , k° to 
find a lower triangular system of linear equations for c\, C2, c n . The lower triangular coefficient 
matrix is invertible (non-singular) because, if this is not the case, a homogeneous solution would 
exist to give a polynomial of a lower degree for the same eigenvalue A. This contradicts to the 
fact that the set {(n + l)(n + 2)} ne N includes only simple eigenvalues. Therefore, a unique value 
for (ci,C2, •••,c n ) exists for a given n. All eigenvectors are linearly independent since polynomials 
of different degrees defined on Z are linearly independent. The set of all eigenvectors gives a basis 
of eigenvectors in the space of analytic functions on Z. 

Finally, we will prove that the basis of eigenvectors for the linear eigenvalue problem (|34p on 
Z m with m ^ 2 is given by {Po,Pi, ...,P m _2j-, which corresponds to the first (m — 1) eigenvalues 
A G {2, 6, ...,m(m — 1)}. This follows from the fact that each polynomial P j is nonz6ro on Irn for 
j G {0, 1, m — 2} in the sense of 

^1^(^)1/0, jG{0,l,...,m-2}. (37) 

By a contradiction, assume that condition (|37p is false, that is Pj(k) has (m — 1) roots on R. 
However, j < (m — 1) and by the Fundamental Theorem of Algebra, Pj(k) = for all k G Z, 
which is a contradiction. Therefore, condition (|37p is satisfied. Furthermore, since polynomials 
{Po, Pi, Pm-2} correspond to distinct eigenvalues, these eigenvectors are linearly indepen- 
dent and form a basis of eigenvectors on X m . This imply that all other polynomials in the set 
{Pj}j^ m _i are linearly dependent from {Po, Pi, P m _2} on I m , which means, in view of dif- 
ferent degrees and distinct eigenvalues, that Pj are identically zero on X m for all j ^ m — 1. 
Therefore, the basis of eigenvectors for the linear eigenvalue problem (|34h on Ifn with vn ^ 2 is 
given by {P ,Pi,...,P m -2}- ■ 

We note that the polynomials Pj(k) in the proof of Lemma [6] are even in k G Z for even j and 
odd in k G Z for odd j. This follows from the parity transformations of operators in (|36[) and the 
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explicit form of the linear eigenvalue problem (|34p . For example, let m = 4 so that X4 = {—1,0, 1} 
and compute eigenvectors and eigenvalues of (|34l) explicitly: 



A = 1 : Cfc = Po(k) = 1, 
A = 3: Ck = Pi{k) = k, 

A = 6 : Cfc = W) = k 



2 3 



For the same case m = 4, ^(A;) = /c(/c 2 — 1) so that Ps(fc) = for all k £ X4. 

We finish this section with the explicit asymptotic approximations for 3-solitons (m = 3). By 
the symmetry of system (|27p with m = 3, we understand that 

xi = X2 = \f2ae~ 1 44> ai = —a, 02 = 0, 03 = a, 

where a is a root of equation 

a - 4V2e-V v ^ ie ~ 1 = 0, 
which is expanded asymptotically as 

o = -^=(-21og(E)-log|log(e)|+21og(2) + o(l)), as e -> 0. (38) 
V2 

Comparison with the asymptotic expansion (|29p shows that log(A _1 l) = or v = 1, which 
means that the asymptotic distribution of frequencies (I32p becomes accurate for m = 3 with O(l) 
replaced by o(l). As a result, we find asymptotic expansions of the two frequencies of out-of-phase 
oscillations near the 3-soliton equilibrium state in the form: 



u? = 2 + (-41og(e) - 2 log I log(e)| + 41og(2)) + o(l), 
= 2 + 3 (-4 log(e) - 2 log I log(e) | + 4 log(2)) + o(l). 



These asymptotic results will be tested numerically in Section [6j 



6 Numerical results 

We now compare the asymptotic results with direct numerical results for the existence and spectral 
stability of 2- and 3-soliton configurations. We identify the relevant branches of stationary solutions 
by solving the ordinary differential equation 

-5«"(0 + ^MO + « 3 (0-m«(0 = o. ( 4 °) 

A fixed point method (Newton-Raphson iteration) is used to solve a discretized boundary- value 
problem, after a centered-difference scheme is applied to the second-order derivatives with a typical 
spacing of A£ = 0.025. The resulting solutions v(£) are obtained starting from the corresponding 
linear eigenfunction (with 2- or 3-nodes at the linear limit) and continuation over the values of 
the chemical potential parameter fx is used in order to extend the branch to the large values of \x. 
Note that the existence and spectral stability of the 1-soliton configuration were examined in our 
earlier work in |12j . 



15 



Once the stationary solution is obtained for each value of /i, we linearize around it, using an 
ansatz of the form: 



r) = v(0 + 5 (a(Oe Ar + 6(£)e^) , (41) 

where 5 denotes a formal (small) parameter. The admissible values of A (eigenvalues) are found 
from the condition that (a, b) G L 2 (M) is a solution of the linear eigenvalue problem 

-*a"(0 + kMO " MO + ^(0(20(0 + &(£)) = *Aa(0, 

+ U 2 H0 - Hb(0 + v 2 (0H0 + 26(C)) = -*A6(0. 1 j 

Using again a discretization of the differential operators on the same grid, we reduce (j42j) to a 
matrix eigenvalue problem which can be solved through standard numerical linear algebra routines. 

Our main results are summarized in Figures [Jj2] for the 2-soliton configuration and Figures [3]|4] 
for the 3-soliton case. 

Fig. [T] compares the numerical result (solid line) for the location of zeros of v(£) to the 
asymptotic expansion (|19jl (dash-dotted line) , where the scaling transformation (j2J) has been taken 
into account to translate the results from e to /i by e = (2/x) -1 . One can see that the asymptotic 
expansion yields a highly accurate approximation of the numerical result. This is also evidenced 
by the right panel of the figure comparing the numerical solution v(£) for fi = 17 (solid line) with 
the variational ansatz (dashed line). 

Fig. [2]shows the smallest eigenvalues of the linear eigenvalue problem (|42p obtained numerically 
(solid line). The resulting eigenvalues can be classified into two types. The first one consists of a 
countable set of pairs of purely imaginary eigenvalues that give frequencies of oscillations of the 
ground state. The main result in Gallo & Pelinovsky [5] states that the frequencies of oscillations 
of the ground state rj e are found in the limit e — ► as follows 



\im uj n (e) = \Jln(n + 1), n ^ 1. 

e^O 

Note that uj\(e) = 2 is preserved for any e > thanks to the symmetry of the Gross-Pitaevskii 
equation with a harmonic potential [12J. Using the scaling transformation ([2]), we conclude that 
these frequencies satisfy the asymptotic limit 

lim Im(A) = ^ ra(r l_ + 5 , n> i. (43) 

The asymptotic limits ([43 p are shown on Fig. [2] by dashed lines. 

The second set of eigenvalues consists of only two pairs of eigenvalues and is associated with 
the relative motions of the dark solitons [17]. One pair of eigenvalues corresponds to in-phase 
oscillations with frequencies Im(A) ~ as fx — > oo (or u ~ y/2 as e — > in notations of the linear 
oscillator equation ([25|) ). The other pair of eigenvalues corresponds to out-of-phase oscillations 
and it is characterized by the asymptotic expansion ([22h . The asymptotic predictions for the 
second set of frequencies are shown by the dash-dotted lines. 

The right panel of Fig. [2] shows the real part of the eigenvalues close to the limit of local 
bifurcation at n = | . The instability, which was studied in [TP] , is caused by the resonance between 
the out-of-phase 2-soliton oscillations and the quadrupolar oscillation mode of the ground state. 
Contrary to what is claimed in numerical work of |19| . we can see from Fig. [2] that the instability 
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Figure 1: Left: the equilibrium position of the two dark solitons versus the chemical potential \i. 
The solid line shows the direct numerical result and the dash-dotted line represents the asymptotic 
approximation (|19f) . Right: the solid line shows the numerical solution v(£) for [i = 17, while the 
dashed line represents the corresponding variational ansatz. 

interval is finite and the 2-soliton excited state may be linearly stable for sufficiently large values 
of the chemical potential fj,. 

We note, however, that the frequency ujq(e) of the out-of-phase oscillations of two dark solitons 
given by the asymptotic expansion (|22p grows as e —* 0. As a result, this frequency will coalesce 
with other frequencies w n (e), n ^ 3 associated with oscillations of the ground state as e — > 0. 
Coalescence with the frequency ujs(e) does not produce an instability, because of the different 
parity of the corresponding eigenfunctions. However, coalescence with the frequency u±{e) will 
produce the instability again and it will happen roughly at e ~ e~ 10 . This value of e is too small to 
be confirmed by our numerical results on Fig. [2J This secondary instability of the 2-soliton excited 
state is anticipated in a tiny interval near e ~ e~ 10 , after which the neutrally stable frequency 
ojo(s) will reappear until further such coalescence occurrences arise with frequencies uje(e), wsfe), 
etc. 

Figures [3] and 0] illustrate similar characteristics but for the 3-soliton state. Once again the 
variational prediction given by the asymptotic expansion (|38|l provides a highly accurate estimate 
of the numerical inter-soliton distance a = — a-i = a% — a\. 

On the other hand, in this case, there exist three frequencies associated with the relative 
motions of three dark solitons, whose values can be seen to be in very good agreement with the 
asymptotic expansion (139 j) . Close to the linear limit n = |, there exists two resonances between 
out-of-phase motion of three dark solitons and the corresponding frequencies of oscillations of the 
ground state. The two resonances induce instabilities of the 3-soliton excited states with two finite 
instability bands. 

The above results provide a relatively complete understanding of the statics and dynamics of 
multi-soliton states within Bose-Einstein condensates at least within the Thomas-Fermi limit of 
large chemical potential. This characterization is especially relevant presently given the recent 
experiments of [171 [T8] enabling the observation and robust time-following for large timescales (of 
the order of hundred milliseconds or more) of such states. However, there would be a multitude 
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Figure 2: Left: solid lines indicate the frequencies of linearization around a 2-soliton solution as 
a function of the chemical potential fi. The dashed lines show the asymptotic limits ([43 1) for the 
frequencies around the ground state. The dash-dotted lines indicate the asymptotic predictions for 
the in-phase (lower frequency) and out-of-phase (higher frequency) oscillations of 2 dark solitons. 
Right: real part of the unstable eigenvalue in a finite instability band near the linear limit of 
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Figure 3: Similar to Fig. [I] but for the 3-soliton case. The left panel again shows the equilibrium 
inter-soliton distance (solid: numerical results; dash-dotted: asymptotic approximation), while 
the right shows the numerical prediction (solid) and variational ansatz (dashed) of the 3-soliton 
state u(£) for /i = 17. 
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Figure 4: Same as Fig. [2j but for the 3-soliton case. The left panel shows the numerical frequencies 
(imaginary parts of the relevant eigenvalues) by solid line, the asymptotic limits for the frequencies 
of the ground state by dashed line, and the frequencies of oscillations of three dark solitons by 
dash-dotted line. The right panel illustrates the real part of the unstable eigenmodes arising close 
to the linear limit \i = ^. 

of directions in which it would be relevant to generalize these results, if possible. On the one 
hand, extending them (analytically) to non-polynomial variants of the Gross-Pitaevskii equation 
accounting for the confinement of the condensate across tranvserse directions would be a challeng- 
ing theoretical task. Another equally interesting direction would involve attempting to generalize 
relevant notions in trying to characterize the dynamics of vortex solitons in higher dimensional 
settings. These directions are presently under consideration and corresponding results will be 
reported in future publications. 
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